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ABSTRACT 

HD 17156b is a newly-found transiting extrasolar giant planet (EGP) that 
orbits its G-type host star in a highly eccentric orbit (e ~0.67) with an orbital 
semi-major axis of 0.16 AU. Its period, 21.2 Earth days, is the longest among 
the known transiting planets. The atmosphere of the planet undergoes a 27-fold 
variation in stellar irradiation during each orbit, making it an interesting sub- 
ject for atmospheric modelling. We have used a three-dimensional model of the 
upper atmosphere and ionosphere for extrasolar gas giants in order to simulate 
the progress of HD 17156b along its eccentric orbit. Here we present the results 
of these simulations and discuss the stability, circulation, and composition in its 
upper atmosphere. Contrary to the well-known transiting planet HD209458b, we 
find that the atmosphere of HD 171 56b is unlikely to escape hydrodynamically at 
any point along the orbit, even if the upper atmosphere is almost entirely com- 
posed of atomic hydrogen and H + , and infrared cooling by H^ ions is negligible. 
The nature of the upper atmosphere is sensitive to to the composition of the 
thermosphere, and in particular to the mixing ratio of H 2 , as the availability of 
H2 regulates radiative cooling. In light of different simulations we make specific 
predictions about the thermosphere-ionosphere system of HD17156b that can 
potentially be verified by observations. 

Subject headings: hydrodynamics — planetary systems — instabilities — plasmas 
— infrared: general 

1. Introduction 



HD17156b i s a tra nsiting, extrasolar giant planet (EGP) that was discovered recently 
by iFischer et al.l (120071 ) as a part of the N2K program, which is a survey of metal-rich stars 
intended to identify short-period planets (IFischer et al.ll2005l ). The transit was first detected 
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by a group of amateur astronomers (IBarbieri et al.ll2007l ). belonging to the Transitsearch.org 
network, which conducts photometric follow-u p observations of kno wn exoplanets that have 
been detected with the radial velocity method (ISeagroves et al.ll2003l ). Following on these ini- 
tial detections, improve d orbital and plane tary parameters for HP 171 56b have recently been 
released by two groups. iGillon et al.l (120071 ) and llrwin et al.l (120081 ). both of whom performed 
independent transit photometry on the system with differe nt ground-ba s ed te lescopes. For 
this work, we have mostly adopted the parameters given in IGillon et al.l (120071 ). 



HP 171 56b orbits a GO star at a distance of ~78 pc from Earth. Based on the observa- 
tions of iFischer et al.l (120071 ). the bolometric luminosity of the star is 2.6 Ls U n and its age 



is ~5.7 Gyr. Its effective temperature, T e ff , = 6079 K, metallic ity [Fe/H] = 0.24, radius 



R* = 1.47 Rsun, and mass, M* = 1.2 M Sun - IFischer et al.l (120071 ) note that the absence of 



Ca II H & K line emission in the spectrum of the star indicates low chromospheric activity. 
They measured Shk = 0.15 and R'hk = -5.04, and derived a rotational period Prot = 12.8 
days for the star. These observations have implications on the likely magnitude of the stellar 
XUV emissions that are an important energy source in the upper atmosphere of the planet. 

HP17156b is in many ways a remarkable planet. The semi-major axis of its orbit is 
a = 0.16 AU, and the orbit is highly eccentric, with an eccentricity of e = 0.67. The orbital 
period of P = 21.2 days is the longest among the known transiting planets. Pue to these 
characteristics, the planet faces a 27-fold stellar flux variation during each orbit, as it moves 
from ~0.27 AU at apastron to 0.052 AU at periastron. In addition to its curious orbit, 
HP 17156b is unique in other ways as well. It is over three times as massive as Jupiter, but 
its radius is comparable to that of Jupiter. The average density of the planet is thus p p = 
3.47 pjup, and this makes it one of the densest gas giants known at present. 

Pue to the high eccentricity of its orbit, and the relatively long semi-major axis, 
HP 17156b is not likely to be rotationally synchronised to its host star. However, the strong 
tidal forces between the planet and the star during; p eriastron are likely to have driven it into 
'pseudo-synchronous' rotation ( Barbieri et al.l 120071 ) . 'Pseudo-synchronous' rotation implies 



synchronisation during periastron, but asynchronous rotation elsewhere along the orbit. The 
character of the orbit and the rotation of the planet aro und its axis make H P17156b an inter- 
esting subject for dynamical studies of its atmosphere. Ilrwin et al. simulated the me- 
teorology in the photosphere of the planet by using the climate model of lLangton fc Laughlin 
(120081 ). and used the results to estimate the phase-dependent 8 pm infrared emissions from 
the planet. In this paper, we present three-dimensional simulations of the upper atmosphere 
and ionosphere of the planet, and investigate the stability of the atmosphere against thermal 
evaporation. 



Absorption by atomic hydrogen in the atmosphere of another close-in EGP, HP209458b, 
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observed in the stellar Lyman a emission line during the transit of the planet across its 
host star, indicates that it is surrounded by an extended atmosphere, composed mainly of 
atomic hydrogen and H + (IVidal-Madjar et al.ll2003l . 120081 ) . The observations indicate that 



the atmosphere is escaping hydrodynamically with a minimum mass loss rate of 10 gs . 
This conclusion is also supported by the tentative detection of oxygen and ionised carbon in 
the upper atmosphere of the pl anet (IVidal-Madiar et al.l 12004 ). the detection of hot atomic 



hydrogen i n the thermosphere (jBallester et al 



Yelldl2004 12009 : iTian et alJl2005 



2007) and numerous mo delling studies (eg. 



Garcia Munozll2007l ; lErkaev et al.ll2007l ) . A recent attempt 



to constrain the mass loss from the planet further yielded a rate of (1.1 ± 0.3) xlO gs 



(jSchneiter et al.ll2007l ). This estimate is based on modelling the Lyman a absorption, taking 
into account the interaction of the exosphere with the incoming stellar wind. 

While it is reasonable to suggest that the atmospheres of other close-in EGPs should 
undergo hydrodynamic escape as well, compelling observational evidence only exists for 
HD209458b. In general, the stability of the atmosphere against hydrodynamic escape de- 
pends on the conditions in the upper atmosphere and ionosphere, which are heated mainly 
by the absorption of stellar XUV radiation. Models of the thermosphere can thus be used to 
assess the overall stability of EGP atmospheres. Using a three-dimensional thermosphere- 
ionosphere model for extrasolar giant planets (EXOTIM), we recently published simulations 
of the up per atmospheres of Jup iter- type EGPs orbiting a solar-type host star at different 
distances (IKoskinen et al.ll2007al ). These simulations indicate that infrared cooling from Hjj~ 
ions may play an import ant role in cooling the thermospheres of EGPs. This Hjj" cooling effect 
was modelled earlier by (jYellell2004l . 120061 ) for HD209458b in a context of a one-dimensional 
non-hydrostatic model. Our calculations extended this effect to three dimensions and a 
wider sample of orbital distances. We also identified a sharp stab ility limit for the atmo - 
sphere of a Jupiter- type exoplanet orbiting the Sun at ~0.15 AU (IKoskinen et al.ll2007bl ). 
Within this limit molecular hydrogen breaks down in the upper atmosphere due to thermal 
dissociation, the atmosphere expands, and begins to escape hydrodynamically. Outside the 
limit, the thermosphere is effectively cooled by H^, and the atmosphere does not undergo 
fast hydrodynamic escape. 

In this context, HD17156b is an interesting subject for similar modelling studies. During 
periastron, the atmosphere of the planet should absorb enough XUV radiation from the star 
to undergo hydrodynamic escape, whereas at apastron distances we would not expect this 
to be the case. Of course, this kind of reasoning is oversimplistic. The planet may not spend 
long enough near periastron for the atmosphere to heat up, or it may not have enough time 
to cool down during apastron. As a result, the atmosphere may be stable throughout the 
whole orbit, or it may escape hydrodynamically at every point. Also, HD17156b is much 
heavier than Jupiter (that we used as a template in our earlier work). Consequently, the 
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scale heights in the atmosphere are shorter, and due to higher gravity, particles need more 
thermal kinetic energy to escape the atmosphere. Higher gravity may thus help to maintain 
stability even at close-in distances during periastron. 

In order to make predictions about the upper atmosphere, ionosphere, and potential 
infrared emissions, we have modified EXOTIM to simulate the progress of HD17156b along 
its eccentric orbit. The details of the new model are discussed in Section [2j Section [3] 
describes the results of the modelling along with some caveats, and Section @] discusses the 
implications of our results in the context of potential observations and relates them to earlier 
work on this subject. 



2. The Model 

2.1. Basic equations 

The model solves the three-dimensional Navier-Stokes equations of continuity, momen- 
tum and energy by explicit time-integration, assuming that the thermosphere is dominated 
by H 2 , He and H and that it is in hydrostatic equilibrium throughout. The basic equations 
are solved on a non-inertial Eulerian, corotating spherical grid, using spherical pressure co- 
or dinates. The details of t he equations, and the model in general, are described extensively 
in iKoskinen et al.l (j2007al ) and are thus not repeated here. However, since the publication 
of that paper, we have added a number of improvements and modifications to the model, 
which are introduc ed below. Some of thes e changes were also described in the Supplemen- 



tary Information of IKoskinen et al.l (l2007bl ). It should be noted that the resulting equations 
resemble the primitive equations of meteorology, with the exception that they include terms 
due to molecular diffusion, which are important in the rare gas of the upper atmosphere. 
Also, our equations fully include the curvature terms arising from spherical coordinates. 

As before, we solve for the horizontal winds from the momentum equations, accounting 
for advection, geopotential gradients, the Coriolis force and viscosity due to molecular diffu- 
sion. Once the horizontal velocities are known, we solve for the vertical velocities from the 
continuity equation. The temperatures are solved from the energy equation, which accounts 
for the advection of energy around the atmosphere by circulation, the absorption of stellar 
X-rays and EUV (XUV) radiation, infrared cooling (by H3 ), heating and cooling due to 
adiabatic contraction and expansion, vertical and horizontal conduction of heat and viscous 
heating. The basic equations of motion are coupled to continuity equations for neutral mass 
fractions that account for vertical and horizontal advection, vertical diffusion and neutral 
chemistry. The numerical solution is based on finite-difference approximations of spatial 
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derivatives and forward-time integration. This scheme is prone to numeri cal instabilities , 
which we filter out by using the two-step smoothing element described by I Shapiro! (Il970l ). 
This element has been adapted for spherical geometry and vector equations, where appro- 
priate. We only apply smoothing in the horizontal direction, and use the filter every 144 s 
of runtime on temperatures and winds, and every 36 s of runtime on mass mixing ratios of 
the neutral species. Generally, the timestep in the simulations varies between 1-2 s. 



2.2. Composition 

As an improvement on the previous model, we have updated the photochemistry and 
calculated new values for the absorption cross sections of H 2 , He and H throughout the 
XUV spectrum. The new reaction list with appropriate references is shown in Table [TJ 
Instead of using coefficients of heat conduction and molecul ar diffusion for atomic hyd rogen 
given by lAchilleos et al.l (119981 ) . we have adopted values from lBanks fc Kockartd (119731 ). The 
coefficients f or other species a re as b efore. We note that the changes do not affect the results 
presented in iKoskinen et al.l (l2007al ) significantly. 



We determine ion densities by calculating photoionisation by stellar XUV radiation, and 
model the subsequent photochemistry by using the reactions listed in Table [T] and assuming 
photochemical equilibrium. In doing so, we effectively ignore ion transport, thus assuming 
that ion lifetimes are short. Also, we calculate the overall neutral density from the ideal gas 
law, assuming that photoionisation has negligible direct impact on neutral mass fractions. 
This assumption is valid as long as ion densities are much smaller than neutral densities. 

In general, the number density of ion species i, n, is given by: 



drii 
~dt 



Pa + LiUi 



where Pi and Lj are the production and loss rates of species i, respectively. For a tidally 
locked planet in photochemical equilibrium, we have Pi + Liiii = 0, and the resulting series 
of non-linear, simultaneous equations ca n be solved by using Newton-Raphson iteration 
(IKoskinen et al.l l2007al ; iPress et al.l Il992l ). However, if the planet rotates asynchronously, 
the rate of change in density rii is not zero in the Eulerian, corotating frame. For such 
planets, our previous model used a crude time-integration technique: 



m(t + 5t) =n(t) + (Pi + Lim)5t (2) 
This method is unstable for short ion lifetimes, which are typical for Hjj", Hj~, and He + . In 
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general, it requires unrealistically short time steps, and it also forced us to assume that all 
is immediately turned into H^. Thus we were unable to model the densities self- 
consistently for asynchronously rotating planets. Clearly, such a limitation would have been 
unacceptable for planet in 'pseudo-synchronous' rotation such as HD17156b. 

Equation ([I]) can be integrated with respect to time by using the method of integrating 
factors, and thus we used an improved time-integration formula: 

m(t + St) = ^(1 - e~ LiSt ) + rn(t)e- LiSt (3) 
Li 

As long as the time step St is kept reasonably short, this method produces accurate ion 
densities for asynhronously rotating atmospheres in photochemical equilibrium. It also allows 
us to fully include H^" in the calculations. 



2.3. XUV Heating 



Previously the model included stellar radiation between 5.0 and 105 nm. We have 
extended this range to cover X-rays, and now include radiation from 0.1 to 105 nm. We note 
that for XUV emissions typical of the current Sun, the inclusion of X-rays in our input has 
almost no impact on the properties of the thermosphere compared to models that do not 
include X-rays. Some authors have proposed that X-rays have a significant impact on the 



atmo spheres of EGPs orbiting young G stars (eg. iCecchi-Pestellini et al.l 120061 ; iPenz et al. 
20081 ). and this may well be the case because a more significant fractio n of the total XUV 
flux from young stars is in X-rays than from mature solar- type stars (IRibas et al.l 120051 ). 



However, it is not true for EGPs orbiting mature main sequence G stars such as the Sun. 

We assume that the XUV flux from HD17156 b is roughly similar to solar flux, and use 
solar maximum fluxes from the SOLAR2000 model (ITobiska et al.ll2000l ) in these simulations. 
This assumption is not likely to be accurate, and considerably variation in the XUV flux 
is possible. Unfortunately, there are no observations or models of the XUV emissions from 
HD17156. The intensity of the emissions is related to th e strength of the stellar magnetic 
dynamo, which depends on the rotation rate of the star (IRibas et al.l 120051 ). The period of 
rotation for HD17156 is half that of the Sun, so the XUV flux could be slightly higher, but 
this is speculative. In the absence of better estimates, we have little choice but to adopt 
solar emissions for the time being. Despite this uncertainty, the models will still help us to 
develop an understanding of the relative dynamics and structure in the thermosphere, and 
perhaps the predictions can be tested with future observations. 
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It should be noted here that we did not perform detailed energetic calculations to work 
out the heating rate due to the absorption of stellar XUV energy. Instead, we assumed 
that 50 % of the absorbed energy is the rmalised in a H?-d ominated atmosphere. This is 
consistent with calculations performed by lWaite et al.l (119831 ) for Jupiter's uppe r atm o spher e 
that account for all other energetic processes apart from Hjj~ cooling. Following lYelld (120041 ). 
we assumed a 10 % heating efficiency in those layers of the thermosphere that are dominated 
by atomic hydrogen. 



2.4. Infrared Cooling 



The only infrared cooling source included in our model is due to emissions. We 
assume that the upper atmosphere is optically thin in the infrared and ignore radiative 
transfer. Thus any emitted infrared radiation escapes directly to space or to the lower 
atmosphere, and in both cases contributes to the cooling of the thermosphere. We neglect 
the possible presence of hydrocarbons or carbon monoxide near the lower boundary of the 
model that may in reality complicate this picture. 

In order to calculate the total emission rate in local thermodyn amic equilibrium (LTE), 
we employ the complete Hjj" line list published by lNeale et al.l (I1995I ). However, the frequency 
of intermolecular collisions in the upper atmospheres of gas giants is generally too low to 
maintain LTE conditions. Previously, we estimated non-LTE emission rates by multiplying 
the LTE rates by a pressure-dependent correction factor, wh ich was derived experim entally 
by fitting to Galileo observations of Jupiter's thermosphere (IKoskinen et al.ll2007al ). As an 
improvement, we now base the correction factor on detailed balance calculations. 

It would be computationally impossible to perform detailed balance calculations for over 
three million line transitions included in the line list, at every grid point in the model, 
and during each time step. In stead, we have sing led out prominent vibrational transitions 
included in an earlier line list of iDinelli et al.l (119921 ) that account for most of the LTE infrared 
emissions. We determined both the LTE and non-LTE emission rates for these transitions, 
adopting detailed balance calculations for the non-LTE emissions, and calculated the ratio 
of the total non-LTE rates to LTE rates for a simple grid of temperatures and densities that 
encompasses the likely values within the model. We then created a table of non-LTE/LTE 
ratios for different temperatures and densities. As the model runs, it interpolates from this 
table bilinearly in order to work out the non-LTE/LTE ratio for each grid point. The final 
emission rate is worked out by multiplying the total LTE emission rate, evaluated by using a 
fourth-order fitted polynomial to the values given by lNeale et al.l (119951 ). by this ratio. This 
method essentially assumes that the non-LTE correction is the same for the complete line list 
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as it is for the reduced line list of iDinelli et al.l (119921 ) . The resulting ratios are listed for the 
substellar point of a typical model run, together with temperatures and press ures, in Table El 



For com parison, experimental correction factors given by equation (15) of iKoskinen et al. 
( I2007al ) are also shown. 



2.5. The Orbit 



Figure[T]is a crude illustration of the orbit of HD 17156b, which also shows the orientation 
of the orbit with respect to an observer on Earth. The transit is offset clo ckwise from 
peria stron by 31.1 degrees. The inclination of the orbit is i = 85.4 degr e es (IGillon et al. 



2007 ). although this figure has bee n revised to 8 6 .5 deg rees by Irwin et al. ( 20081 ). It is not 
clear if a secondar y eclipse exists. IGillon et al.l (120071 ) suggest that it should be a partial 
grazing eclipse, but llrwin et al.l (120081 ) ascribe only a 9.2 % chance for this to be case. The 
probability of a full secondary eclipse is even lower at 6.9 %. In case secondary eclipse does 
occur, the angular separation between the antitransit and periastron is 148.9 degrees. 

Simulating the planet's progress on its eccentric orbit is done by varying the irradiation 
conditions continuously as the planet moves along the orbit. For this purpose, we need to 
work out the true anomaly, which is the angular distance of the planet from periastron, and 
the distance of the planet from the star as a function of time. Also, as we assume that the 
planet is locked in 'pseudo- synchronous' rotation around its axis, we need to account for the 
planet's spin by varying the position of the star in the corotating frame of the planet. 



The 'pseudo-synchronous' spin angular velocity is given by (jHutlll98ll ): 



n 



sp 



1 + (15/2)e 2 + (45/8)e 4 + (5/16)e e 



n 



(4) 



3e 2 + (3/8)e 4 ](l-e 2 ) 3 / 2 

where e is the eccentricity, and Q is the mean orbital angular velocity. Using e ~0.67 yields 
Q sp ~5.6 Q. This means that during one full orbit the planet spins around its axis 4.6 times 
in the star's frame of reference. 



In terms of the orbital angular velocity during periastron, Q p , the 'pseudo-synchronous' 
spin angular velocity can be expressed as: 



1 + (15/2)e 2 + (45/8)e 4 + (5/16)e e 



[l + 3e 2 + (3/8)e 4 ](l + e) 2 " v (5) 

Again using e ~0.67 yields Q sp ~0.818 Q p . In other words, during periastron the planet's 
spin is slower than the orbital angular velocity, and while passing the periastron the planet 
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should therefore revolve 'backwards' with respect to the star compared to normal, faster spin 
at other parts of the orbit. This behaviour is a curious feature of 'pseudo-synchronisation' 
and it causes a peculiar jitter in the position of the star in the planet's sky near periastron. 

The orbital mean anomaly is given by: 



M(t) = T (t-T) (6) 

where P is the orbital period and T is the time of periastron. This is simply the angular 
distance traversed by the planet in time (t — T) if the orbit was circular. Mean anomaly 
can be related to the eccentric anomaly by using Kepler's second law, which states that the 
radius vector from the host star to the planet sweeps out equal areas in equal time, and the 
geometric properties of an ellipse. The relation between the mean and eccentric anomalies 
is given by Kepler's equation: 



M = E - esin(E) (7) 

where E is the eccentric anomaly. By using the general equation of an ellipse and geometric 
relations, it is straightforward to show that the true anomaly, 9, is related to the eccentric 
anomaly by: 

cos(E) — e 

cos(6) = y ' (8) 
1 — ecos{E) 

Once true anomaly is known, it is easy to solve for the orbital distance as a function of time 
by using: 

d(t) = - Q(1 ~ & 2 „ (9) 
w 1 + ecos[6(t)] v 1 

where a is the semi-major axis of the orbit. 

Kepler's equation <£T§ is a transcendental equation, which cannot be solved analytically. 
Fortunately, it can be solved easily by using Newton-Raphson iteration. For this purpose, 
the equation is first written as: 



f(E) = E - esinE - M = 



If the initial guess at the solution is given by E g , then the progressive corrections to this 
solution are given by: 
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E - E9 ~1W) 

where f'{E g ) is the differential of / with respect to E. The iteration proceeds until an 
acceptable level of convergence is achieved. In general, the solution to Kepler's equation 
converges very quickly and only a few iterations are required. Analytically, it can be shown 
that for e <0.99 and an initial guess of ir, convergence is guaranteed. 

The change in the local hour angle in time St is given by the numerical difference of the 
spin angular velocity and orbital angular velocity multiplied by St: 



SH a = [Q sp - Q(9)]St (10) 

where Q(9) is the orbital angular velocity, which depends on the true anomaly of the planet's 
position. Note that near periastron, where Q(9) is faster than Q sp , SH a can become negative. 

The above suite of equations allows us to model the stellar irradiation on the atmosphere 
at every point along the orbit. During every time step the model calculates a new value for 
the mean anomaly. This value is then converted into true anomaly by using equations (j7j) 
and (jHl). The distance to the host star is then calculated from equation (Q and this distance 
is used to obtain the dilution factor for the stellar XUV flux. Equation (TTOT) is used to 
calculate the position of the star in the planet's sky. Given that this procedure is repeated 
every timestep, it proceeds surprisingly swiftly and does not add much to overall computation 
time. 

Figure [2] shows the orbital true anomaly versus time for the orbit of HD17156b, Figured 
shows the orbital distance versus true anomaly and Figure 0] shows how the hour angle 
develops during one orbit. All simulations start from apastron, with 9 = 180°, and the 
local zenith for the hour angle calculation is set initially at the substellar point. The plots 
illustrate that the orbital solution is working as it should and that it makes physical sense. 
Also, they demonstrate the fact that the orbital angular velocity is faster during periastron, 
implying that the planet spends most of the time completing the 'far-side' of the orbit. 



2.6. Boundary Conditions and Planetary Parameters 



We used planetary and orbital parameters that are listed in Table [3] together with 
the other input pa rameters for the simulations. We used parameters originally released by 
Gillon et al.l (120071 ) in their preprint. These parameters have been subsequently revised, and 
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they also differ from those released by llrwin et al.l (120081 ). However, given all the uncertainties 



in the model, and the uncertainties in the published parameters, the differences are not likely 
to be particularly important. The numerical solution is evaluated on a grid of 36 evenly 
spaced longitude points, 31 evenly spaced latitude points, and 28 pressure levels ranging 
from 2 /ibar to 0.04 nbar with a uniform spacing of 0.4 scale heights. The calculations do 
not extend to the poles and instead, the field variable values are interpolated over to the 
polar latitude circles (31,30,1,2). We assume that the thickness of the lower atmosphere is 
entirely negligible compared to the radius of the planet, and thus place the zero-altitude 
reference point at our lower boundary at 2 /ibar. 

We assume that the lower boundary temperature is constant at T =520 K. This is 
close to the equilibrium temperature of a gas giant with a Bond albedo equal to 0.3 under 
solar irradiation near the apastron of HD17156b's orbit. Using a constant lower boundary 
temperature for a planet like HD17156b may not be realistic as this temperature is likely 
to vary along the orbit. Unfortunately it is not clear how the temperature varies at the 2 
/ibar level as the irradiation level changes. In order to avoid speculation we have chosen to 
use a constant lower boundary temperature for our reference models. However, as the effect 
on our model could be significant, we also produced simulations where we varied the lower 
boundary temperature. The results are discussed in Section l3~5l 



The initial composition comes from a model of Jupiter's auroral region (IGrodent et al. 



20011 ) that effectively fixes the mixing ratios of the neutral species at the lower boundary of 
the model. The lower boundary composition for a planet like HD17156b is highly uncertain. 
It is not clear if the dominant carbon-bearing molecule in the lower atmosphere should be 
CO or CH 4 , and there are no constraints on the mixing ratios of heavier molecules in general. 
Thus we have chosen not to speculate on hydrocarbon or oxygen photochemistry, and instead 
assumed that the mixing ratios of heavier molecules in the thermosphere are negligible. This 
assumption is not as crude as it may initially appear, because due to molecular diffusion the 
mixing ratios of molecules heavier than He should be very small inside our pressure range. 

The mixing ratios of H and He near our lower boundary are also uncertain, which is why 
our first instinct was to use jovian values (~2 x 10 -4 for H, and 0.056 for He). However, these 
values may well be inappropriate for a planet that is so different from Jupiter as HD17156b. 
For instance, atomic hydrogen in the upper atmosphere is created as a by-product of the 
complex oxygen and hydrocarbon photochemistry, which is powered by stellar UV radiation, 
and is expected to be taking place in t he stratosphere-mesosphere region of close-in gas 



giants. Calculations of iLiang et al.l (120031 ). based on such photochemistry, demonstrate that 
the mixing ratio of atomic hydrogen near the /ibar level could be around 1 % for a close-in 
giant like HD209458b (with a = 0.045 AU). They indicate that for a variety of chemistries, 
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the production of atomic hydrogen is mainly limited by the strength of the UV flux, which 
varies 27-fold for HD17156b. Based on these arguments, we expect the mixing ratio of atomic 
hydrogen to be less than 1 % for HD17156b at the /xbar level because most of the time the 
planet is farther from the host star than HD209458b, and thus it should receive less UV 
radiation. 

While our model is not sensitive to small variations in the mixing ratio of He at the 
lower boundary, the amount of atomic hydrogen could be crucial because it regulates the 
availability of H 2 and thus limits the formation of hht. In order to examine the effect of 
changing the lower boundary mixing ratio of atomic hydrogen, we have generated both 
Jupiter-type models and models where the mixing ratio of atomic hydrogen was raised to 1 
% at the lower boundary at the expense of the mixing ratio of H 2 . 

For all of our simulatio ns, we assume zero win ds at the lower boundary This point 



was extensively discussed in iKoskinen et al.l (j2007al ) , where we also demonstrated that the 



inclusion of winds and temperature variations at the lower boundary does not affect the 
conditions in the upper atmosphere significantly. This insensitivity is due to the very strong 
forcing by the stellar XUV radiation in the middle-to upper thermosphere. As an upper 
boundary condition, we assume that vertical gradients reduce to zero in the two uppermost 
layers, vertical winds vanish at the upper boundary, and that the neutral species are in 
diffusive equilibrium in the two uppermost levels. 



3. Results 

3.1. Evaporation of the Atmosphere 

The stability of the atmosphere against thermally driven hydrodynamic escape is deter- 
mined at the exobase, which is defined as the altitude where an escaping particle undergoes 
approximately one collision within one atmosph eric scal e heig ht. A handy measure of sta- 



bility is the thermal escape parameter, given by iHuntenl (119731 ): 



A = GM p m = ,V esc. 2 

where m is the mass of the escaping constituent, k is Boltzmann's constant, T c is the tem- 
perature at the critical level (exobase), and z c is the altitude of the critical level. This is 
basically the square of the ratio of the escape velocity to the most probable thermal velocity, 
u, at the exobase. Once the thermal escape parameter becomes less than about 1.5 for H or 
less than about 3.0 for H 2 , the thermal kinetic energy at the critical level becomes greater 
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than the gravitational potential energy and the atmosphere begins to escape hydro dynami- 
cally. If A >3, Jeans escape dominates. Of course, the transition between Jeans escape and 
hydrodynamic escape is not abrupt, and some overlap is to be expected for low values of A. 
However, at A >30, even Jeans escape is fairly insignificant. 

The escape velocity from HD17156b is 108 kms -1 , which compares to about 60 kms -1 
for Jupiter. Thus the temperature in the upper atmosphere of HD17156b must be signifi- 
cantly higher for hydrodynamic escape to take place than what is required for a Jupiter-like 
planet. Figure shows the globally averaged temperatures for HD17156b at the upper 
boundary of the model versus orbital true anomaly (starting from -180° at apastron). The 
results are shown for two types of simulations, one with lower boundary mixing ratios from 



Grodent et al.l (120011 ) (Exo-1) and one with the mixing ratio of atomic hydrogen fixed at 1 



at the lower boundary (Exo-2). 



For Exo-1, the exobase temperature is around 2,000 K at apastron, and it rises to over 
3,000 K at periastron. The temperature peaks around periastron and reaches a minimum 
near apastron. The temperature difference at apastron between successive orbits is ~0 K, 
indicating that the model repeats the same behaviour during each orbit and is therefore 
in a kind of steady state. Remarkably, this steady state establishes itself after only three 
simulated orbits. The upper thermosphere is dominated by molecular hydrogen, and the 
thermal escape parameter varies between 690 during apastron and 410 during periastron, 
rendering the atmosphere stable throughout the whole orbit and even Jeans escape negligible. 
As expected, the thermosphere is relatively thin, with the upper boundary altitude ranging 
from 840 km during apastron to 1,100 km during periastron. 

For Exo-2, the exobase temperatures are considerably higher, over 10,000 K throughout 
the orbit. The minimum temperature is around 15,000 K, reached at 9 = -145°, and the 
maximum temperature is 34,000 K, reached after the periastron at 9 = 76°. Atomic hydrogen 
and H + are virtually the only species in the upper thermosphere, and between the minimum 
and maximum temperature regions the thermal escape parameter ranges from 40 to 15, 
respectively. This implies that near periastron, Jeans escape could be significant and even 
some bulk flows are possible. However, according to the criteria based on the thermal escape 
parameter, the atmosphere should not undergo fast hydrodynamic escape at any point along 
the orbit. The upper boundary altitude ranges from 13,100 km (1.2 R p ) near the temperature 
minimum to 27,300 km (1.4 R p ) near the temperature maximum. The mass loss rate from 
Exo-2, based on thermal Jeans escape, peaks around 10 6 gs _1 at around 9 =76°. However, 
this mass loss rate is not sustained throughout the orbit, and instead it drops to less than 
10~ 2 gs -1 around apastron. 

The temperature variation between successive apastron passes for Exo-2 is ~130 K after 
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fifteen simulated orbits. This means that the simulations have not reached exact steady 
state, alhough it is reasonable to assume that they are near steady state. 130 K is not 
a particularly significant figure compared to the temperature of the upper boundary (over 
15,000 K), and the difference gets smaller during each orbit. Computational constraints 
mean that we cannot run the Exo-2 simulation to exact steady state within a reasonable 
time period, and thus it is possible that the model keeps heating up slowly during each orbit 
until it reaches conditions that would lead to hydrodynamic escape. Given the rather swift 
passage through the periastron, however, we think that this is unlikely. 

The drastic difference in the thermal characteristics of the two simulations of HD17156b, 
Exo-1 and Exo-2, arises from the differences in the composition of the thermosphere. Fig- 
ure M shows the substellar mixing ratios of atomic hydrogen during apastron for the two 
simulations. For Exo-1, the mixing ratio of atomic hydrogen is relatively low at all levels, 
rising to about 6 % near the exobase. At the upper boundary, the maximum mixing ratio 
is about 10 % near the dawn terminator whereas the minimum is about 5.5 % near the 
dusk terminator. For Exo-2, atomic hydrogen becomes the dominant species between 1 and 
0.1 /ibar (altitude of 50-300 km). Above that level, atomic hydrogen and H + are virtually 
the only species in the thermosphere. This is because, initially, the higher mixing ratio of 
atomic hydrogen at the lower boundary in the Exo-2 simulation leads to a reduction in the 
atomic hydrogen concentration of the thermosphere through molecular diffusion. This limits 
the formation of fht, which reduces the significance of radiative cooling compared to the 
Exo-1 simulation. As a result, the temperatures become high enough to dissociate molecular 
hydrogen thermally, and this process quickly converts most of the upper atmosphere into 
atomic hydrogen that is subsequently ionised by the intense XUV radiation. It should be 
noted, though, that for HD17156b the transition region between the H2 and H atmospheres is 
located deeper in the thermosphere than it would be for a Jupiter-type planet or HD209458b 
with similar lower boundary conditions because of the higher gravity on HD17156b. 

Figure [7] contrasts the apastron H3 densities and infrared cooling rates at the substellar 
points of Exo-1 and Exo-2. The substellar column densities of Hjj~ in the two models are 
3.2 xlO 16 m~ 2 and 5.0 xlO 14 m~ 2 , respectively, and in Exo-2 hLt virtually disappears at 
pressures lower than 0.6 /ibar (above 150 km). This difference in densities makes itself 
evident in the infrared cooling rates, that are prominent for Exo-1 but negligible for Exo- 
2. Contrasting the total infrared cooling rates to the total XUV heating rate allows us to 
estimate global cooling functions for the models. Depending on the orbital positions, the 
cooling function for Exo-1 is 72-91 % while for Exo-2 it is merely 0.1-0.24 %. The cooling 
function thus depends on the composition and photochemistry of the models, and the degree 
to which H_2 is dissociated thermally. 
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It is unfortunate that the conditions in the thermosphere depend on the lower boundary 
composition, as the mixing ratios there are difficult to constrain. On the other hand, this 
feature opens up interesting new avenues for observations. For instance, an observation 
of H Lyman a absorption during transit could be used to constrain the properties of the 
upper atmosphere of HD1 7156b in the same way as the detection of the extended hydrogen 
cloud around HD209458b (IVidal-Madjar et al.l 120031 ) has allowed for the characterisation of 
that planet's atmosphere. We predict that for the upper atmosphere of HD17156b to be 
dominated by atomic hydrogen, it should extend to more than 1.4 R p during transit near 
periastron. If an extended hydrogen atmosphere is not detected, then infrared cooling must 
be taking place and the thermosphere is dominated by H 2 . Thus, even if B.^ emissions are 
too faint to be detected directly, we could infer their existence indirectly and also constrain 
the mixing ratios of H and H2 in the lower atmosphere. 



3.2. Circulation 

Figure [8] shows the temperature and circulation at the upper boundary of Exo-1 at four 
different orbital positions: during apastron, a quarter orbit after apastron (6 = -153°), during 
periastron, and a quarter orbit after periastron {6 = 153°). If secondary eclipse takes place, 
the fourth orbital position we have considered here is near (although not exactly at) the 
longitude of the antitransit, wheras the periastron models are a fairly good representation 
of the transit conditions. 

At all four orbital positions, the dayside is clearly warmer than the night side. At 
apastron, the dayside temperature is 2,200-2,300 K, and the temperature peaks at the equator 
near dawn. The night side temperature, with a minimum close to the dawn terminator, is 
1,650-1,750 K. The diurnal temperature difference is thus 500-600 K. The horizontal winds 
originate in the dayside, and blow accross the terminator to the night side. In the night 
side, the Coriolis force and the geopotential gradients drive the eastward wind towards the 
equator, into a stream that flows across the antistellar point and faces the westward wind 
from the dayside near the dawn terminator. Also in the night side, there are two high- 
latitude vortices that direct the easterly wind from the dayside into the westerly equatorial 
jet. The maximum zonal wind speed is 1.1-1.3 kms -1 at high latitudes. Convection is 
directed upwards in the dayside, accompanied with adiabatic expansion of the atmosphere, 
and downwards in the night side, accompanied with adiabatic contraction of the atmosphere. 
The maximum vertical flow speed is only a few ms" 1 . The thermal escape parameter for the 
average atmosphere at the upper boundary is 689, indicating negligible thermal evaporation. 
An escaping particle undergoes approximately 5 collisions within one scale height above the 
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upper boundary, implying that the 0.04 nbar level is slightly below the exobase. 

At 9 = -153 degrees the distance of the planet from the star is 0.22 AU. The upper 
boundary temperature and wind pattern is qualitatively identical to the apastron model. 
However, with added heating, the dayside temperatures are higher, at 2,300-2,450 K. The 
night side temperature is 1,650-1,850 K, implying a diurnal temperature difference of 600-700 
K. The maximum zonal wind speed is also slightly higher at 1.3-1.6 kms -1 . The thermal 
escape parameter is 659, indicating that thermal evaporation is insignificant. 

At periastron, where the orbital distance is 0.052 AU, the near-synchronisation of the 
planet's spin leads to a more pronounced diurnal temperature difference. Between 9 = -90° 
and 9 = 90°, only one hemisphere of the planet is exposed to stellar irradiation. The dayside 
temperature during periastron is around 4,350 K, while the night side temperature is around 
2,200 K. Thus the diurnal temperature difference is over 2,000 K. Large-scale circulation 
is qualitatively similar to other orbital positions, but the winds are faster, with maximum 
speeds of 2.6-2.8 kms" 1 . Vertical advection is also more rigorous, and the maximum vertical 
wind speed is 8-15 ms" 1 . However, despite the added heating, the thermal escape parameter 
is 411, and evaporation is still slow. The enhanced XUV heating is balanced by effective 
Hjj" cooling, and before the thermosphere has time to heat up significantly, the planet moves 
away from periastron. At 9 = 153°, the upper boundary is almost identical to the 9 = -153° 
model. There are some very slight differences between these models, but those are mostly 
limited to the lower thermosphere. 

Deeper in the thermosphere, the situation is different. Figure [9] shows the tempera- 
ture and circulation near 55 nbar (altitude of 250-300 km) for Exo-1 at the same orbital 
positions as above. In all cases, the temperature is nearly uniform. With the exception of 
periastron, the circulation is characterised by a broad eastward, circumplanetary jet. During 
apastron, the temperature varies between 1,530 and 1,580 K, and the equator is everywhere 
slightly warmer than its surroundings, with temperature peaking in the night side. The wind 
speed in the equatorial jet ranges from 160 to over 220 ms^ 1 . Cyclonic polar vortices circle 
around the minimum temperature regions near the poles. In general, the circulation in the 
lower thermosphere is characteristic of Coriolis-driven dynamics, that smoothes out diurnal 
temperature variations. 

At 9 = -153°, the temperature varies between 1,540 and 1,610 K. The equator is still 
warmer than its surroundings, and there are two temperature peaks, one near the substellar 
point and one in the night side. The temperature minima are again located near the poles 
near the dawn terminator, and the cyclonic polar vortices are centred on those. In the 
equatorial jet, the wind speed ranges from 100 to 220 ms -1 . 
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During periastron the horizontal temperature variations are more pronounced, with the 
temperature ranging from 1,840 to 2,020 K, but the details of the temperature distribution 
are confusing. The warm substellar region is surrounded by a cooler ring along the termi- 
nator, but the night side is almost as warm as the dayside. The winds diverge from the 
substellar point, blowing towards the night side. An eastward jet develops in the night side, 
but it encounters the easterly wind at dawn, and cannot encircle the planet at this pressure 
level. The two high-latitude vortices direct easterly flows into the equatorial jet. The maxi- 
mum wind speed is around 800 ms" 1 . The circulation is qualitatively similar to the topside 
circulation, and this must be a consequence of the near-synchronisation during periastron. 

At 9 = 153°, the temperature ranges from 1,580 to 1,660 K. Overall, the lower thermo- 
sphere is slightly warmer than it was at the symmetric position at 9 — -153°, and it is in the 
process of cooling down after the periastron passage. The circulation is again characterised 
by the circumplanetary jet and the polar vortices. The wind speed in the equatorial jet 
ranges from 100 to 400 ms" 1 . 

For comparison, Figure [TD] shows the temperature and circulation at the upper boundary 
of Exo-2 for the same orbital positions. The horizontal uniformity of the temperature in these 
simulations is remarkable, and the wind speeds are very slow. During apastron, the dayside 
temperature is 17,620-17,630 K, while the night side temperature is around 17,530 K. Thus 
the diurnal temperature difference is only of the order of 100 K. The winds blow from the 
dayside to the night side, and in the night side some of the eastward flow is directed into an 
equatorial jet that faces the easterly wind from the dayside near dawn. High-latitude vortices 
turn some of the easterly flow around, and direct it into the equatorial jet. In the dayside, 
the temperature peak is shifted slightly eastward from the substellar point by circulation 
and rotation. The winds are slow, with a maximum speed around 200 ms -1 . An escaping 
particle suffers ~14 collisions within one scale height at the upper boundary, indicating that 
the 0.04 nbar level is located below the exobase. However, the P-T profile is isothermal at 
pressures lower than ~0.7 nbar, and thus the temperature at the critical level should be the 
same as at our upper boundary. The thermal escape parameter at this boundary is 33, and 
the mass loss rate based on Jeans escape is ~10~ 2 gs _1 . 

At 9 = -153°, the model is generally cooler than at apastron. The dayside temperature 
is around 15,180 K, while the night side temperature is around 15,030 K. The fact that the 
model cools down despite increasing XUV fluxes implies that a slow cooling process takes 
place after periastron that continues even after apastron. The circulation is qualitatively 
identical to the apastron model, but the maximum wind speed is slightly higher at 340 
ms -1 . The thermal escape parameter increases to ~39 at the upper boundary, and the mass 
loss rate is ~10~ 4 gs~ x . 
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Moving from 9 = -153 to periastron the model is heated by the steadily increasing XUV 
flux. Between these two positions, the flux is multiplied by a factor of 18. During periastron, 
the dayside temperature is 29,200-29,500 K, with a maximum on the equator, shifted slightly 
eastward from the substellar point. The night side temperature is around 26,800 K, implying 
that the diurnal temperature difference near transit should be of the order of 2,500 K. This 
implies that during periastron the XUV flux drives the model, instead of dynamics and slow 
cooling. The winds blow from the dayside to the night side, converging near the antistellar 
point, with a maximum speed of 2.8-3.2 kms -1 . The thermal escape parameter is 19, which 
indicates that the atmosphere is still in the thermal Jeans escape regime. The mass loss rate 
is ~10 4 gs -1 . 

By the time the planet reaches 9 = 153 degrees, the thermosphere has started cool- 
ing down. At this point, the dayside temperature is around 23,430 K and the night side 
temperature is around 23,270 K. Thus the diurnal temperature difference is of the order of 
150 K. As the upper thermosphere slowly cools, energy is advected around the planet and 
this mixing smoothes out the diurnal temperature difference that develops during periastron. 
Curiously, the day-night circulation is still symmetric about an axis connecting the substellar 
point to the antistellar point. The winds have slowed down to less than 200 ms" 1 , and the 
thermal escape parameter has increased to ~23 at the upper boundary. The mass loss rate 
is ~2 xlO 3 gs" 1 . 



3.3. Thermal Structure and Energy Balance 

Figure [IT] shows the substellar P-T profiles for the Exo-1 series at the four orbital 
positions discussed in the previous section. Outside periastron the P-T profiles are isothermal 
at pressures lower than about 0.2 nbar. Generally the temperature increases with altitude, 
and the gradient is steepest in the lower thermosphere at pressures higher than 0.3 /ibar. 
Overall, the temperatures are highest during periastron, and lowest during apastron. At 
9 = -153° and at 9 = 153° the P-T profiles are identical in the upper thermosphere, but 
the lower thermosphere is warmer after periastron. During periastron the temperature rises 
sharply with altitude towards the upper boundary due to added XUV heating. 

Figure [T2l shows the substellar heating and cooling rates by volume for the Exo-1 simu- 
lations at the same orbital positions. Generally, in the middle and upper thermosphere (at 
p < 0.1 /ibar), the XUV heating is effectively balanced by H^" cooling. However, at pressures 
lower than 1.0 nbar, this cooling is less important due to non-LTE conditions, and the heating 
is balanced mainly by conduction of heat downward. In the lower thermosphere, the heating 
is balanced by the 'triad' of vertical advection, adiabatic expansion, and conduction. Due 
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to relatively low temperature, Hjj" cooling is not that significant near the lower boundary. In 
the night side (not shown), the thermosphere is mainly heated by vertical downwelling and 
adiabatic contraction. The heating is balanced by vertical conduction. 

In Figure [T2| the solid line shows the net heating rate (the sum of all the energy equation 
terms) at each orbital position. The net heating rate is close to zero in the upper thermo- 
sphere, while variations take place in the lower thermosphere. The radiative timescale in the 
thermosphere is relatively short, and in those regions where H^" cooling is significant, any 
added heating is quickly balanced by an adjustment in the cooling rate. This explains the 
rapid convergence of the Exo-1 simulations to steady state after only a few orbits. Circula- 
tion in the upper thermosphere is driven by uneven stellar irradiation, which explains the 
relatively steep diurnal temperature gradient. 

In the lower thermosphere, on the other hand, net heating occurs during periastron. By 
the time the planet reaches 9 = 153° (0.22 AU), the lower thermosphere has started to cool 
down, due to conduction, advection and adiabatic expansion. The cooling continues until the 
model reaches 9 = -153°. After this the model starts to heat up again as it moves towards 
periastron. The lower part of the model, instead of being radiation- driven, is dynamics- 
driven, and thus its response to the varying XUV fluxes is not as rapid as that of the upper 
thermosphere. Also, the diurnal temperature differences tend to disappear at high pressures, 
and as we have seen, the circumplanetary wind appears. 

We define the radiative cooling function as the ratio of the total infrared cooling rate to 
the total XUV heating rate. This cooling function is relatively high for the Exo-1 simulations. 
During apastron, it is 78 %, and the total Hjj~ infrared emission rate (in all spectral lines) is 
6.1 xlO 14 W. At 9 = -153°, the cooling function is 72% and the emission rate is 8.5 xlO 14 
W. During periastron, the corresponding figures are 85 % and 1.7 xlO 16 W, respectively. 
Near the antitransit, and thus the possible secondary eclipse, the cooling function is 91 % 
and the emission rat e is 1.1 xlO 15 W. As a proxy for potentially detectable spectral lines, 
Shkolnik et al.1 (120061 ) estimated that the total Hjj" emission rate should be of the order of 10 17 - 



10 19 W if emissions are to be detected with current ground-based technology. This estimate 
is based on a sample of F, G, K, and M stars within a distance of ~40 pc. Unfortunately, 
even for these relatively high cooling functions, the predicted emission rates fall well short 
of the detection limits. 

Figure [13] shows the substellar P-T profiles for the Exo-2 simulation at the same four 
orbital positions. Outside the periastron, the P-T profiles are isothermal at pressures lower 
than about 0.7 nbar, and during periastron the profile is isothermal at pressures lower than 
about 0.1 nbar. In the lower thermosphere, the temperature increases steadily with altitude. 
The temperature in the upper thermosphere increases as the planet moves towards periastron 
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from 9 = -153°. The heating goes on for a while after periastron, but by the time the planet 
reaches 9 = 153°, the outer layers have started to cool down. Curiously, this does not apply 
to the region between 3.0 and 100 nbar, where the model is actually warmer at 9 = 153° 
than it is during periastron. Towards apastron, the whole thermosphere cools down and this 
cooling continues until the planet reaches 9 = -153° again. 

Figure [14] shows the corresponding volume heating and cooling rates for the above P-T 
profiles. The XUV heating (by volume) is concentrated near the bottom of the thermo- 
sphere. This is because, compared to H2, atomic hydrogen is not particularly effective in 
absorbing XUV radiation. Also, the heating efficiency in the outer envelope is only 10 %. 
The radiation thus penetrates to the H 2 -dominated bottom layers, where it is also absorbed 
by He, and where the heating efficiency is higher. However, the temperature is higher in 
the upper thermosphere, because the lower thermosphere is much denser and heats up slug- 
gishly whereas even relatively inefficient heating is enough to produce high temperatures in 
the outer layers. 

Generally, the stellar heating is balanced by vertical conduction, vertical advection and 
adiabatic expansion. In the night side, again, the heating is by downwelling and adiabatic 
contraction. Hjj" cooling is negligible compared to other cooling mechanisms in the dayside, 
and generally it only occurs at the bottom of the thermosphere. As indicated by the net 
heating rates, the thermosphere heats up as it approaches periastron after 9 = -153°, and the 
heating first occurs in the lower thermosphere. The model reaches temperature maximum 
after periastron, and by the time the planet reaches 9 = 153°, the lower and upper ther- 
mospheres have started cooling down. Around the 100 nbar level XUV heating and vertical 
conduction are still heating the thermosphere at 9 — 153°, and this is reflected in the P-T 
profiles above. 

The radiative cooling functions are low for the whole orbit. At apastron, the cooling 
function is a meagre 0.23 % and the total Hjj" emission rate is 4.9 x 10 11 W. This is three 
orders of magnitude less than the emission rate for the Exo-1 apastron model. The cooling 
functions for 9 = -153°, periastron, and 9 = 153° are 0.24, 0.097, and 0.19 %, respectively. 
The corresponding Hg~ emission rates are 7.6 xlO 11 W, 5.6 xlO 12 W, and 6.3 xlO 11 W. 

3.4. Composition 

Obviously, there are huge differences in thermospheric composition between the two 
scenarios, Exo-1 and Exo-2, with potentially observable consequences. Figure [T5l shows the 
substellar number densities of the neutral species during apastron for both models. For 
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Exo-1, H2 dominates at all levels. Atomic hydrogen overtakes helium around 0.1 /ibar as the 
second most prominent neutral species. Qualitatively, this picture prevails throughout the 
whole orbit. 

For Exo-2, atomic hydrogen is the dominant species at pressures lower than 1 /ibar, and 
at higher altitudes the concentrations of the other species are negligible. Intriguingly, helium 
dominates over H 2 at pressures between ~4 and 600 nbar. This is because the decline of H 2 
in the upper thermosphere is not merely due to molecular diffusion, but it is significantly 
enhanced by thermal dissociation. 

Figure [TBI shows the substellar electron density profiles for Exo-1 and Exo-2, both dur- 
ing periastron and during apastron, and Figure [17] shows the corresponding ions density 
profiles for H + , , and hLt. For Exo-1, the electron density during apastron increases with 
altitude in the lower thermosphere, and reaches a peak value of ~1.8 xlO 13 m -3 around 
36 nbar. In general, the electron density is close to the density of H + , which is by far the 
dominant ion in the thermosphere. The electron density drops steeply in the night side in 
the lower thermosphere, but in the outer layers, the densities are fairly uniform around the 
whole atmosphere. This is because the lifetime of H + against recombination in the upper 
thermosphere is of the order of 40 hours, and the timescale for planetary rotation is similar. 
The relatively long lifetime of H + indicates that the assumption of photochemical equilib- 
rium is inaccurate in the upper thermosphere, where ion transport is likely to be important. 
However, the lifetimes of other ions are relatively short, from seconds to minutes, and even 
the lifetime of H + is shorter in the lower thermosphere. The situation in this respect is thus 
similar to Jupiter, where H + densities also deviate from photochemical equilibrium in the 
upper thermosphere. The H3 profile in Figure [T7J has two distinct peaks, one between 0.7 
and 5.0 nbar, and the other near the lower boundary. The substellar column densities of H + 
and are 7.9 xlO 18 m~ 2 and 3.2 xlO 16 m~ 2 , respectively. 

During periastron, there is a strong peak in the electron density profile near the 0.6 
/ibar level, with a density of 8.9 xlO 13 m -3 . This peak coincides with a H + peak in the 
lower thermosphere. Intriguingly, the electron and H + densities near the upper boundary 
are actually lower during periastron than they are during apastron. There are two main 
reasons for this. First, the temperature in the upper boundary region is higher during 
periastron, and this leads to a lower overall density. Second, the near-synchronisation of 
the planet's spin leads to enhanced day-night circulation that increases the mixing ratio 
of n 2 in the dayside. The horizontal winds are fuelled by upwelling in the dayside and 
downwelling in the night side. Upwelling brings up gas from the lower thermosphere, where 
the mixing ratio of H2 is higher, and thus replenishes H2 concentrations at higher altitudes. 
Conversely, in the nightside the mixing ratio of atomic hydrogen is embellished as H2 sinks 
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with downwelling. This kind of circulation is an important factor amplifying the cooling 
function in the thermospheres of tidally locked planets. During apastron, the mixing of the 
atmosphere is more effective, as the relative rotation of the planet with respect to the star is 
faster. This leads to more uniform mixing ratios of atomic hydrogen, which do not exhibit 
notable diurnal variations. For the periastron model, the substellar column densities of H + 
and Hjj~ are 2.0 xlO 19 m~ 2 and 6.3 xlO 17 m~ 2 , respectively. The hht densities are quite a 
bit higher than at apastron, and the mixing ratio of ion is also higher in middle and upper 
thermosphere. 

For Exo-2, during apastron, there are two peaks in the electron density profile at the 
bottom of the thermosphere and around 36 nbar, both with the roughly equal density of 6.3 
xlO 13 m~ 3 . Overall, the electron densities are higher than in the Exo-1 simulation. As H + is 
virtually the only ion in the thermosphere, the electron density profile is effectively identical 
to the H + profile. The lifetime of H + in the upper thermosphere is around 90 hours, and 
a significant plasma density persists in the night side. This is such a large deviation from 
photochemical equilibrium that ion transport cannot be ignored. Unfortunately this means 
that our density plots are not likely to be accurate. 

At pressures lower than about 0.6 //bar, there is virtually no Hjj" in the Exo-2 simulation. 
Very small amounts of are present, but the thermosphere does not contain enough H2 
to convert this into Hj~. In fact, the density of He + is generally higher in the thermosphere 
than the density of Kj\ The substellar column densities of H + and are 7.9 xlO 20 m -2 
and 5.0 xlO 14 m~ 2 , respectively. 

During periastron, the peak electron densities are much higher, around 4.5 xlO 14 m~ 3 . 
The relative ion concentrations are similar to the apastron model. At pressures lower than 
~0.7 nbar, H + is actually the dominant species. In fact, the mixing ratio of H near the upper 
boundary, in relation to H + , is only 30 %. We have not made any provision for this in our 
simulations, which assume that the neutral density can be calculated from the ideal gas law, 
and that ionisation does not affect neutral mass fractions directly. This is a further reason 
to treat the Exo-2 ion densities with suspicion, at least during periastron. Elsewhere along 
the orbit, the density of H + tends to be much smaller compared to H. 

One last uncertainty related to the calculations of neutral and ion densities in the Exo-2 
series is the fact that the numerical solution for the neutral mass fractions uses H 2 as the 
'dominant' species, although the equations themselves do not assume a dominant species and 
diffusion is calculated for three major species instead of minor species diffusion. However, 
the numerical solution evaluates the mass fractions for He and H, and then deduces those 
from unity in order to obtain the mass fraction for H2. This approach leads to numerical 
errors in the H-dominated part of the thermosphere. The problem is limited to the relatively 
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thin transition region between the H2-dominated layers anf the H-dominated layers. Above 
this region the mixing ratio of atomic hydrogen is unity for all practical purposes. 

3.5. Variable temperatures 

In order to explore the effect of variable lower boundary temperatures on our simulations, 
we generated a model, labelled Exo-3, which is otherwise identical to Exo-1 but where 
we varied the lower boundary temperature assuming that it is given by the equilibrium 
temperature of the planet (with a Bond albedo of 0.3) at different locations along the orbit. 
We updated this temperature at every time step during the simulation. Figure [TS] shows 
the lower boundary temperature as a function of true anomaly and contrasts this with 
the globally averaged temperature at the upper boundary of the Exo-3 simulation. The 
equilibrium temperature varies from ~490 K at apastron to just over 1,100 K at periastron. 
This kind of variation, which is basically symmetric around the periastron, is unlikely to be 
realistic as the response of the lower atmosphere to varying stellar irradiation is likely to be 
more complex. 

However, as Figure [TBI indicates, varying the lower boundary temperature has almost no 
effect on the upper boundary temperatures. As it turns out, changing the lower boundary 
temperature does not affect the qualitative nature of the circulation and composition of the 
Exo-1 simulation either. Figure [18] also shows the substellar P-T profiles from the two simu- 
lations during periastron. At pressures between about 0.6 /zbar and 1 nbar the P-T profiles 
coincide. Above the 1 nbar level, the Exo-3 model is slightly warmer, but the difference in 
temperature is barely noticeable. In the lower thermosphere the temperatures are obviously 
different due to the different lower boundary conditions. The XUV heating rates and the IR 
cooling rates in the two models are almost identical. In the lower thermosphere the dayside 
cooling effect due to vertical advection is slightly reduced in the Exo-3 model, while the 
cooling effect due to adiabatic expansion is enhanced. Also, due to the shallow vertical P-T 
profile near the lower boundary, the vertical conduction term produces slight heating near 
the lower boundary of the Exo-3 simulation. In conclusion, it seems that the effect of the 
higher lower boundary temperature in the Exo-3 simulation is dwarfed by the fact that the 
temperatures in the thermosphere are largely determined by radiative balance, and most of 
the XUV radiation is absorbed above the lower boundary. 

We did not run a separate simulation to explore how the varying lower booundary 
temperature would affect the Exo-2 simulations, but as the upper boundary temperature in 
those simulations is consistently above 10,000 K, it is unlikely that a variation of a few 100 
K at the lower boundary of the model would have a significant effect on the results. 
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4. Conclusions and Discussion 

We have generated three-dimensional simulations of the upper atmosphere of HD17156b, 
by using an improved vers ion of the coupled ther mosphere-ionosphere model for extrasolar 
giant planets (EXOTIM) jKoskinen et all l2Q0Tal fbh . By adopting the model for this EGP 



that orbits its G-type host star on a highly eccentric orbit (e =0.67), we have been able 
to make predictions about the nature and stability of its upper atmosphere, thermospheric 
dynamics and ion concentrations. These predictions have observable consequences and if 
followed up, they could prove very useful in probing the conditions in the atmosphere of this 
exotic world. 

We have found that the composition of the thermosphere is sensitive to the composition 
of the underlying layers, and in particular to the mixing ratio of atomic hydrogen in the 
lower atmosphere. Due to molecular diffusion, we do not anticipate a significant presence 
of heavy molecules s uch as CO, CH4 and N 2 , that would lead to extensive and complicated 



photochemistry (eg. iGarcia Munozl 120071 ). in the upper atmosphere. The availability of H 2 
regulates the Hjj~ cooling function in the thermosphere, and thus different mixing ratios of 
atomic hydrogen can lead to very different conditions in the upper atmosphere. Unfortu- 
nately, the mixing ratio of atomic hydrogen in the lower atmosphere is poorly constrained. 
Due to this uncertainty, we have generated two different types of HD17156b simulations: 
one with lower boundary mixing ratio of atomic hydrogen appropriate for Jupiter, and one 
with atomic hydrogen dominating the upper atmosphere. 

In both cases, we find that the atmosphere of HD17156b remains stable everywhere 
along the orbit and does not begin to escape hydrodynamically. This is despite the fact 
that a planet with a circular orbit at the periastron distance of HD17156b (0.052 AU) from 
the host star would almost certainly develop a planetary wind similar to that observed on 
HD209458b. Evaporation is not as significant on HD17156b because the planet does not 
spend long enough within 0.1 AU from the parent star for hydrodynamic escape to begin. 
Also, it is three times as massive as Jupiter, and thus the escape velocity and the thermal 
energy required for escaping the atmosphere is relatively high. 

Different thermospheric compositions, and thus different radiative cooling functions, 
lead very different conditions in the upper atmospheres. If the thermosphere is dominated 
mostly by H 2 , the H3" cooling function varies between 70 and 90 % along the orbit, and 
the total infrared emissions are of the order of 10 15 -10 16 W. In this case the extent of the 
atmosphere is negligible, and even the upper atmosphere stays relatively cool (2,000-3,000 
K). The thermal escape parameter at the upper boundary ranges from 400 to 600, and 
this renders any kind of evaporation negligible. If, on the other hand, the thermosphere is 
dominated by H and H + , the radiative cooling function is only 0.1-0.24 %, and the total 
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infrared emissions are of the order of 10 12 -10 13 W. In this case the upper atmosphere is hot 
(15,000-30,000 K) and during periastron, the atmosphere expands beyond 1.4 R p . Depending 
on the orbital position, the thermal escape parameter ranges from 14 to 40, indicating that 
the atmosphere is in the Jeans escape regime. 

These predictions can potentially be verified by observations. The infrared emissions 
from Hjj~ are potentially observable in the future, and they would provide the most direct 
clue to the properties of the upper atmosphere and its interaction with the host star. Un- 
fortunately, the predicted inte nsities are so low, tha t detection cannot be made with current 



ground-based technology (eg. IShkolnik et al.l 120061 ) . The extended atomic hydrogen cloud 



could possibly be observed by monitoring H Lyman a absorption during transit, which oc- 
curs near periastron. It is not clear, however, if such an observation is possible in the near 
future because the instrument that performed a similar measurement on HD209458b, i.e. 
STIS onboard HST, is no longer available. 

Orbital dynamics and the planet's spin affect three-dimensional simulations in many 
important ways. Circulation in the upper thermosphere is characterised by day-night winds, 
driven by a persistent diurnal temperature contrast. This behaviour is enhanced during peri- 
astron, when the planet is in near-synchronisation with the star. In the lower thermosphere, 
temperatures are fairly uniform horizontally and the circulation is characterised by an east- 
ward, circumplanetary jet. In the H-dominated simulations, temperatures are generally more 
uniform and the wind speeds are slow. 

As thermal escape is not as significant for HD17156b as it is for, say, HD2009458b, 
one wonders if non-thermal escape processes or tidal forces might come into play instead. 
Non-thermal escape arises from charge exchange with the impinging stellar wind, photodis- 
sociation or dissociative recombination, acceleration of ions along open magnetic field lines, 
sputtering of neutral particles or ion pick-up by the stellar wind, and impact ionisation and 
dissociation by stellar wind particles. Quantifying these effects for a planet like HD17156b 
is difficult, because there is no information about its magnetic field or the properties of the 
stellar wind in its vicinity. Esti mates of non-therm al mass loss for a planet like HD209458b 
are of the order of 10 8 -10 9 gs -1 ( lErkaev et al.ll2005l ). If similar values are to be expected for 
HD 17156b, then evaporation from the planet is mostly non-thermal. 



Adopting the formula presented by lErkaev et al.l (120071 ). the approximate Roche lobe 
distance during periastron fro m HD17156b is 10-11 R p . This m eans that the occurence of 
so-called geometrical blow-off (ILecavelier des Etangs et al.l 120041 ). driven by the tidal forces 
between the planet and the star, is very unlikely. In fact, it may be safe to ignore tidal forces 
in estimating the mass loss for HD17156b. 
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Table 1. Chemical reactions used by the model. 





Reaction 




Ks.te a ' h 


Reference 


1. 


H 2 + hu -: 


• H J+ c 






Yan et al. ( 1998) 


2. 


H 2 + hu -* 


• H++H + c 






Chung et al. fl993),Duiardin et al. (1987) 


3. 


R 2 + hu -: 


• H++H++ e 






Duiardin et al. f 1987) 


4. 


H +hu -* 


H++ c 






Hummer et al. (1963) 


5. 


He +hv -> 


He++ e 






Yan et al. f 1998) 


6. 


H 2 + M -> 


2H + M 


1.5 x 10" 


9 exp(-4.Se4:/T) 


Baulch et al. (1992) 


7. 


2H + M - 


. H 2 + M 


8.0 x 10" 


33 (300/T)°.6 


Ham et al. (1970) 


8. 


H++ H 2 (> 


> 4) -> H++ H 


1.0 x 10" 


■ 9 exp(-2.19e4/T) 


Yelle (2004). estimated 


9. 


H++ H 2 - 


♦ H++ H 


2.0 x 10" 


9 


Thread & Huntress (1974) 


10 


. H++ H 2 + M -> H++ M 


3.2 x 10" 


29 


Kim fc Fox (1994) 


11 


. He++ H 2 


-♦ H++ H + He 


1.0 x 10" 


9 exp(-5700/T) 


Moses fc Bass (2000) 


12 


. He++ H 2 


H++ He 


9.35 x 10 


-15 


Anicich (1993) 


13 


. H++ H - 


+ H++ H 2 


2.0 x 10" 


9 


Yelle (2004). estimated 


14 


. H++ H - 


-> H++ H 2 


6.4 x 10" 


10 


Kapras et al. ( 1979) 


15 


. H++e -» 


H +hu 


4.0 x 10" 


12 (300/T e ) - 64 


Storey & Hummer (1995) 


16 


. H++e - 


H + H 


2.3 x 10" 


8 (300/T e ) 04 


Auerbach et al. ( 1977) 


17 


. He++e - 


► He +hv 


4.6 x 10" 


12 (300/T e ) - 64 


Storev & Hummer (1995) 


18 


. H++e -c 


H 2 + H 


2.9 x 10" 


8 (300/T e ) - 64 


Sundstrom et al. ( 1994) 


19 


. H 3 ++e -> 


H + H + H 


8.6 x 10" 


8 (300/T e ) - 64 


Datz et al. (1995) 



a Photo-ionisation rates are calculated explicitly by using the photoionisation cross sections given in the 
references. 

b Two-body rates are given in cm 3 s — 1 and three-body rates are given in cm 6 s _1 . The electron temperatures 
are assumed to be the same as neutral temperatures. 
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Table 2. Non-LTE correction factors 



Pressure (nbar) 


Temperature (K) 


Detailed balance 


Experimental 


9nnn 
zuuu 




i nn 


1 nn 

l.UU 




ooo 


i nn 


1 nn 

l.UU 


oyy 


1119 
1 1 1Z 


i nn 

l.UU 


1 nn 

l.UU 


602 


1382 


0.998 


1.00 


404 


1614 


0.986 


1.00 


271 


1732 


0.971 


1.00 


181 


1792 


0.956 


1.00 


121 


1849 


0.934 


1.00 


82 


1925 


0.907 


1.00 


55 


2008 


0.870 


0.878 


36 


2083 


0.824 


0.782 


25 


2147 


0.767 


0.696 


16 


2201 


0.702 


0.618 


11 


2252 


0.630 


0.548 


7.4 


2311 


0.556 


0.484 


5.0 


2382 


0.480 


0.426 


3.3 


2465 


0.404 


0.374 


2.2 


2555 


0.330 


0.327 


1.5 


2654 


0.258 


0.284 


1.0 


2766 


0.194 


0.245 


0.67 


2899 


0.138 


0.210 


0.45 


3074 


0.0916 


0.179 


0.30 


3314 


0.0554 


0.150 


0.20 


3615 


0.0303 


0.124 


0.14 


3919 


0.0155 


0.101 


0.09 


4144 


0.0082 


0.0794 


0.06 


4254 


0.0050 


0.0602 


0.04 


4254 


0.0036 


0.0428 
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Table 3. Planetary and orbital parameters for HD 171 56b 



Planetary parameters" 




M p {M Jup ) 


3.1 


Ftp (Rjup) 


0.96 


Orbital parameters" 




Eccentricity 


0.67 


Semi-major axis (AU) 


0.16 


Inclination (°) 


85.4 


Period (days) 


21.2 


Longitude of periastron (°) 


121 


Lower boundary conditions 




To (K) 


520.0 


Po (Pa) 


0.2 


Gravity (ms — 2 ) 


87.0 


Mixing ratio of H (Exo-1) 


2 xl0~ 4 


Mixing ratio of H (Exo-2) 


0.01 



a With the exception of the planetary ra- 
dius, the orb ital and plane t ary pa rameters were 
taken from IGillon et al.l | |2007| ). Note that 
iGillon et al.|[|2007l) have recently published and 
estimated radius of 1.23 ^n'* 7 Rjun for t he 
planet. The value given bv llrwin" et all ||2008| ) is 
1.01 ±0.09 Rj up . Thus, given the error ranges, 
these two estimates coincide. We have used a 
value that a grees with the estimates given by 
Ilrwin et alj feOOSl) . 
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Fig. 1. — The orbit of HD1 7156b. Note that the orbit is not drawn to scale and is meant 
for illustration purposes only. The planet is shown at different positions near periastron to 
illustrate 'pseudo-synchronous' rotation. During the periastron the same side of the planet 
always faces the star, but as the planet moves farther from true anomaly of 90°, it begins to 
rotate asynchronously. The angle uo marks the longitude of the periastron. 
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Time (doys) 



Fig. 2. — Orbital true anomaly versus time (in Earth days). The orbital period of HD17156b 
is approximately 21.2 Earth days. Periastron is located at 9 — 0. The planet spends only 
~2 days within 0.1 AU, passing through periastron. 
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Fig. 3. — Orbital distance (in AU) versus true anomaly. The orbital period of HD17156b is 
approximately 21.2 Earth days. Periastron is located at 6 — 0. 
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Fig. 4. — Hour angle versus orbital true anomaly. Here the hour angle traces the angular 
separation of a fixed point on the planet's surface and the substellar point. Zero and 360 
degrees correspond to the substellar point. 
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Fig. 5. — Globally averaged upper boundary (at 0.04 nbar) temperatures for two different 
simulations with different lower boundary mixing ratios of atomic hydrogen and radiative 
cooling functions. The solid line shows temperatures for a model with a cooling function 
of 70-90 % and the dotted line shows temperatures for a model with a cooling function of 
0.01-0.24 % (see text). 
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Fig. 6. — The mixing ratio of atomic hydrogen at the substellar point of Exo-1 (solid line) 
and Exo-2 (dotted line) during apastron. 



-38- 



(o) 




-10 -5 5 10 15 -3x10"'° -2x10"'° -1x10"'° 

Log[Number density (m"*)] Infrored cooling (Wm"*) 



Fig. 7. — H3 densities and the corresponding volume infrared cooling rates for (a) Exo-1 
and (b) Exo-2 (see text) during apastron at the substellar point. 
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Fig. 11. — Substellar P-T profiles for the Exo-1 simulations at apastron (solid line), 9 = -153° 
(dotted line), periastron (dashed line), and 6 = 153° (dash-dotted line). 
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Fig. 12. — Substellar energy equation terms (per unit volume) for the Exo-1 simulations at 
(a) apastron, (b) 6 = -153°, (c) periastron, and (d) 6 = 153°. 
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Fig. 13. — Substellar P-T profiles for the Exo-2 simulations at apastron (solid line), 9 = -153° 
(dotted line), periastron (dashed line), and 9 = 153° (dash-dotted line). 
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Fig. 14. — Substellar energy equation terms (per unit volume) for the Exo-2 simulations at 
(a) apastron, (b) 6 = -153°, (c) periastron, and (d) 6 = 153°. 





Fig. 15. — Substellar neutral number densities during apastron for (a) Exo-1, and (b) Exo-2. 
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Fig. 16. — Substellar electron density profiles for (a) Exo-1 and (b) Exo-2, during apastron 
(solid line) and periastron (dotted line). 
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Fig. 17. — Substellar ion density profiles for Exo-1 during (a) apastron and (b) periastron, 
and for Exo-2 also during (c) apastron and (d) periastron. 
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Fig. 18. — (a) Globally averaged upper boundary (at 0.04 nbar) temperature (solid line) 
for the simulation with a variable lower boundary temperature (Exo-3), and the equilibrium 
temperature of a gas giant with a Bond albedo of 0.3 (dotted line), (b) Substellar P-T 
profiles for the Exo-1 (solid line) and Exo-3 (dotted line) simulations at periastron. 



